##
# Name:
#    GWA_analysis_functions.R
# Author:
#    Floriane Simonet
# Mail:
#    simonet.floriane@gmail.com
##

# necessite les packages : car et gap

################# VUE D'ENSEMBLE DES DONNEES et elimination de SNPs #####################

#  Creation d'une liste de SNPs sur les chromosomes X, Y et MT et/ou comportant les combinaisons alléliques A/T, T/A, G/C et C/G
del_ATGC_chrXYMT <- function(inbfile=NULL,del_AT_GC=T,del_chrXYMT=T,outputfileATGC="AT-GC.txt",outpufilechrXYMT="chrXYMT.txt"){
	if(!is.null(inbfile)){

		bim <- read.table(file=paste(inbfile,".bim",sep=""))

		# Creation d'une liste de SNPs sur les chromosomes X, Y et MT
		if(del_AT_GC==T){
			chrXYMT <- bim[bim$V1 > 22,"V2"] 
			write.table(as.data.frame(chrXYMT),file=outpufilechrXYMT,quote=F,row.names=F) 
			cat(paste("see heterosomes'list : ",outpufilechrXYMT,"\n",sep=""))
		}

		# suppression A/T, T/A, G/C et C/G
		if(del_chrXYMT==T){
			bim$bin <-ifelse((bim$V5=="A" & bim$V6=="T"),0,1) # delete A/T
			bim$bin <-ifelse((bim$V5=="T" & bim$V6=="A"),0,bim$bin) # delete T/A
			bim$bin <-ifelse((bim$V5=="C" & bim$V6=="G"),0,bim$bin) # delete C/G
			bim$bin <-ifelse((bim$V5=="G" & bim$V6=="C"),0,bim$bin) # delete G/C

			delATCG <- as.character(bim[bim$bin==0,"V2"])
			write.table(unique(delATCG),file=outputfileATGC,quote=F,row.names=F)
			cat(paste("see ATGC SNPs : ",outputfileATGC,"\n",sep=""))
		}
	}
}


#####


#creation d'une liste de SNPs à eliminer en fonction de hardy-weinberg, "minor allele frequency" et du pourcentage de missing
snp_selection <- function(inbfile=NULL,HW=NULL, miss=NULL, MAF=NULL,outputfile="SNPtoDel.txt"){
	if(!is.null(inbfile)){

		data <- NULL
		tmp_list <- NULL
		delHW <- NULL
		delmiss <- NULL
		delMAF <- NULL

		# creation d'un nombre aléatoire servant de prefixe
		randnb <- round(runif(10,min=0,max=1000))
		prefix <- randnb[1]

		cat("create SNPs list to be deleted : \n")

		if(!is.null(HW)){
			# creation d’un fichier de desequilibre de Hardy Weinberg (hwe.log, hwe.hh, hwe.hwe)
			system(paste("plink --noweb --bfile ",inbfile," --hardy --out ",prefix,"hwe",sep=""))
			hwe <- read.table(file=paste(prefix,"hwe.hwe",sep=""),h=T)
			
			# recherche du status des individus du fichier(cas, contrôle ou cas/contrôle)
			fam <- read.table(file=paste(inbfile,".fam",sep=""))
			fam2 <- fam[fam$V6 != -9,]
			study <- length(unique(fam2$V6)) # study = 1 (cases or controls), study = 2 (cases and controls) 

			if(study==1){
				if(unique(fam2$V6) == 1){hwe <- hwe[hwe$TEST=="UNAFF",];}
				if(unique(fam2$V6) == 2){hwe <- hwe[hwe$TEST=="AFF",];}

				del_hwe <- hwe[hwe$P < HW & !is.na(hwe$P),] 

				# creation de la liste de SNPs a supprimer
				delHW <- c(as.character(del_hwe$SNP))
			}
			else if(study == 2){
				aff <- hwe[hwe$TEST=="AFF",]
				unaff <- hwe[hwe$TEST=="UNAFF",]
				tmp <- merge(aff,unaff,by=c("CHR","SNP"))

				tmp_aff <- tmp[tmp$P.x < HW & !is.na(tmp$P.x),"SNP"] 
				tmp_unaff <- tmp[tmp$P.y < HW & !is.na(tmp$P.y),"SNP"]

				# creation de la liste de SNPs a supprimer
				delHW <- c(as.character(tmp_aff),as.character(tmp_unaff))
			}
			
			cat("add SNPs (p-value hwe < ",HW,")\n",sep="")
		}

		if(!is.null(miss)){
			# creation d'un fichier de difference de call rate (% de genotypage) entre cas et contrôles (dmiss.log, dmiss.hh, dmiss.missing)
			system(paste("plink --noweb --bfile ",inbfile," --test-missing --out ",prefix,"dmiss",sep=""))	
			dmiss <- read.table(file=paste(prefix,"dmiss.missing",sep=""),h=T)

			if(study==1){
				if(unique(fam2$V6) == 1){del_dmiss <- dmiss[dmiss$F_MISS_U > miss,"SNP"];}
				if(unique(fam2$V6) == 2){del_dmiss <- dmiss[dmiss$F_MISS_A > miss,"SNP"];}

				#  creation d'une liste de  SNPs avec un % de missing  > miss
				delmiss <- as.character(del_dmiss)
			}
			
			if(study==2){
				tmp_aff2 <- dmiss[dmiss$F_MISS_A > miss & !is.na(dmiss$F_MISS_A),"SNP"]
				tmp_unaff2 <- dmiss[dmiss$F_MISS_U > miss & !is.na(dmiss$F_MISS_U),"SNP"]

				# creation d'une liste de SNPs avec un % de missing  > miss
				delmiss <- c(as.character(tmp_aff2),as.character(tmp_unaff2))
			}

			cat("add SNPs (missing < ",miss,")\n",sep="")
		}

		if(!is.null(MAF)){
			# creation d'un fichier de frequence allelique (maf.log, maf.hh, maf.frq)
			system(paste("plink --noweb --bfile ",inbfile," --freq --out ",prefix,"maf",sep=""))
			maf <-  read.table(file=paste(prefix,"maf.frq",sep=""),h=T)

			# creation d'une liste de SNPs consideres comme "non frequent" (< MAF)
			delMAF <- as.character(maf[maf$MAF < MAF,"SNP"]) 
			cat("add SNPs (minor allele frequency < ",MAF,")\n",sep="")
		}


		# concaténation des données et suppression des doublons
		delSNP <- c(as.character(delHW),as.character(delmiss),as.character(delMAF))
		final_list <- unique(delSNP)

		# enregistrement des donnees
		write.table(final_list,file=outputfile,quote=F, row.names=F,col.names=F)
		cat(paste("see results : ", outputfile,"\n",sep=""))

	}
	else { cat("missing inbfile parameter.\n")}
}


#####


# elimination des SNPs en fonction des F_U et F_A (fréquences allélique chez les cas et les contrôles)
del_freq <- function(inbfile=NULL, th_freq=NULL, outputfile="SNPoutfreq.txt"){
	if(!is.null(inbfile) & !is.null(th_freq)){
		system(paste("plink --noweb --bfile ",inbfile," --assoc",sep=""))
		assoc <- read.table(file="plink.assoc",h=T)
		assoc2 <- assoc[na.omit(assoc$F_A) < th_freq,]
		assoc3 <- assoc[na.omit(assoc$F_U) < th_freq,]
		assoc4 <- assoc[na.omit(assoc$F_A) > (1 - th_freq),]
		assoc5 <- assoc[na.omit(assoc$F_U) > (1 - th_freq),]
		del <- c(as.vector(assoc2$SNP),as.vector(assoc3$SNP),as.vector(assoc4$SNP),as.vector(assoc5$SNP))
		write.table(del,file=outputfile,quote=F,row.names=F)
		cat("see results : ",outputfile,"\n")
	}
}

#####

# visualisation des individus en fonction de leur heterozygotie 
heterozygosity <- function(inbfile=NULL,hist_breaks=NULL,outputfile="heterozygosity"){

	diff.het <- function(data){
		res <- as.numeric(data[c(4)]) - as.numeric(data[c(3)])
	}

	if(is.null(hist_breaks)){
		rm(hist_breaks)	
	}

	# creation d'un nombre aléatoire servant de prefixe
	randnb <- round(runif(10,min=0,max=1000))
	prefix <- randnb[1]

	system(paste("plink --noweb --bfile ",inbfile," --het --out ",prefix,"het",sep=""))
	het <- read.table(file=paste(prefix,"het.het",sep=""),h=T)
	het$diffEO <- apply(X=het,MARGIN=1,FUN=diff.het)
	
	write.table(het,file=paste("summary-het",outputfile,".txt",sep=""),quote=F,row.names=F)	

	png(paste("plots-",outputfile,".png",sep=""),height=580,width=920)
	par(mfrow=c(1,2))
	plot(het$diffEO)
	hist(het$diffEO,breaks=hist_breaks)
	dev.off()

	cat(paste("\n see plots-",outputfile,".png and summary-het",outputfile,".txt \n",sep=""))
}

########################## STRATIFICATION ############################

# fonction permetttant de faire le plot de tous les PI HAT, le plot des moyennes de PI HAT par individu et l'histogramme de la distribution des moyennes de PIHAT 
PIHAT.plots <- function(genomefile="plink.genome", famfile=NULL, general_plot_title="PI-HAT", mean_plot_title="average PI-HAT for each individual",distrib_title="distribution",distrib_breaks=100,outputfile="plots_PI-HAT"){

	if(!is.null(famfile)){

		# chargement des donnees
		genome <- read.table(file=genomefile,h=T)
		fam <- read.table(file=famfile)
		genome <- genome[genome$IID1 %in% fam$V2,]
		genome <- genome[genome$IID2 %in% fam$V2,]

		# creation des graphiques concernant les donnees PIHAT
		filename <- paste(outputfile,"png",sep=".")
		png(filename, width=860, height=860)
		par(mfrow=c(2,2))

		# visualisation du nuage de points du PI-HAT des individus
		cat(paste("editing general PIHAT plot with",genomefile," \n",sep=" "))
		plot(genome$PI_HAT, xlab="individuals pair", ylab="PI-HAT", main=general_plot_title)

		# PI_HAT moyen pour chaque individu
		cat("writing averagePIHAT.txt\n")
		mean_PIHAT <- tapply(X=genome$PI_HAT,INDEX=list(genome$IID1),FUN=mean)
		IID <- as.character(names(mean_PIHAT))
		table <- cbind(IID,mean_PIHAT)
		write.table(table, file="averagePIHAT.txt",quote=F,row.names=F)

		# plot des PI_HAT moyen pour rechercher les individus à exclure
		cat("editing means PIHAT plot\n")
		plot(mean_PIHAT, main= mean_plot_title, xlab="individuals", ylab="PI-HAT")
	
		# distribution des donnees
		cat("editing means PIHAT distribution\n")
		hist(mean_PIHAT, main=distrib_title, breaks=distrib_breaks, xlab="mean PIHAT")

		dev.off()
	
		cat(paste("see results :",filename,"\n",sep=" "))
	}
	else{
		cat("\nERROR : missing famfile parameter \n")
	}

}

#####

# creation d'une liste d'individu a exclure de l'analyse
PIHAT.exclusionlist <- function(genomefile="plink.genome", famfile=NULL, th_indiv=NULL, th_mean=NULL, controls=NULL){

	if(!is.null(famfile)){

		# chargement des donnees
		means <- read.table(file="averagePIHAT.txt",h=T)

		# apres visualisation graphique on decide de supprimer les individus qui ont une moyenne de PI_HAT > th_mean
		if(!is.null(th_mean)){
			indiv_out <- means[means$mean_PIHAT > th_mean,]
			# recuperation des FID
			ind <- read.table(file=famfile,h=F)
			del_list <- unique(ind[ind$V2 %in% indiv_out$IID, c("V1","V2")])		

			# creation du fichier d'individus a exclure
			write.table(as.data.frame(del_list), file="indiv_meanpi-hat_to_remove.txt", row.names=F,col.names=F,quote=F)
			cat("see results list : indiv_meanpi-hat_to_remove.txt\n")
		}
	
		if(!is.null(th_indiv)){
			if(!is.null(th_mean)){
				indiv <- means[means$mean_PIHAT < th_mean,]
			}else{
				indiv <- means
			}
			genome <- read.table(file=genomefile,h=T)
			genome2 <- genome[genome$IID1 %in% indiv$IID,]
			genome3 <- genome2[genome2$IID2 %in% indiv$IID,]
			genome4 <- genome3[genome3$PI_HAT > th_indiv, c("FID1","IID1","FID2","IID2")]
			genome4$bin <- ifelse(genome4$IID2 %in% genome4$IID1,1,0) # on ne prend qu'un IID sur 2 
			IIDlist <- genome4[genome4$bin==0, c("FID1","IID1")]
			if (!is.null(controls)){
				addcontrols <- unique(ind[ind$V2 %in% controls,c("V1","V2")]) # ajout de tous les controles
				colnames(addcontrols) <- c("FID1","IID1")
				final_list <- unique(rbind(IIDlist,addcontrols))
				write.table(final_list,file=paste("list_indiv_pihat_sup",th_indiv,".txt",sep=""),quote=F,row.names=F,col.names=F)
			}else{
				write.table(IIDlist,file=paste("list_indiv_pihat_sup",th_indiv,".txt",sep=""),quote=F,row.names=F,col.names=F)
			}
			cat("see results list :",paste("list_indiv_pihat_sup",th_indiv,".txt",sep=""),"\n")
		}
	}
}

#####

# creation d'un fichier phenotype pour les composantes
create_strat_phenofile <- function(mdsfile="strat.mds", outputfile="pheno_strat.txt"){

	pheno <- read.table(file=mdsfile,h=T)
	pheno2 <-pheno[,-c(3)]
	write.table(pheno2, file=outputfile ,row.names=F,quote=F)
	cat("see results : ",outputfile,"\n")
}

#####

# creation de plots montrant l'association des individus aux composantes
plot.indiv.components <- function(nb_comp=NULL, famfile=NULL, phenofile="pheno_strat.txt" ,outputfile="qqplot_components", color=NULL){# color=gender for coloring by gender, color=status for coloring cases/controls

	if(!is.null(nb_comp) & !is.null(famfile)){
		# chargement des donnees
		fam <- read.table(file=famfile)
		strat <- read.table(file=phenofile, h=T)
		col_names <- colnames(strat[-c(1,2)])
		strat <- merge(strat, fam, by.x="IID", by.y="V2")
		strat <- strat[-c(1,2)] 
		colnames(strat) <- NULL
	
		nb_plot <- ((nb_comp)*(nb_comp -1))/2
		nb_page <- ifelse(nb_plot%%4 == 0, nb_plot%/%4, (nb_plot%/%4)+1)
		nb_output <- 1 # compte le nombre de fichiers en sortie
		png_bool <- 0 # permet la gestion du nombre de fichiers png (en fonction du nombre de plots dejà mis sur le png)
		last_plot <-0

		for(indA in 1:(nb_comp -1) ){
			for(indB in (indA+1):nb_comp){
				compA <- strat[c(indA)]
				compB <- strat[c(indB)]
				if(color=="gender") {
					col <- strat[c(length(strat)-1)]
					# creation des parametres de la legende 
					col_legend<- cbind(c(1,2),c("male","female"))
				}
				if(color=="status"){
					col <- strat[c(length(strat))]
					col_legend<- cbind(c(1,2),c("controls","cases"))
				}	
				# creation de la matrice de donnees servant au plot
				tmp <- cbind(compA,compB,col)	
				if(nb_page > 1){
					if(png_bool == 0){png(paste(outputfile,nb_output,".png",sep=""), width=860, height=860); par(mfrow=c(2,2))}#on ouvre un nouveau png tous les 4 plots
					# creation du plot
					plot(tmp$compA, tmp$compB, main="plot components", xlab=col_names[indA], ylab=col_names[indB],col=tmp$col)
					legend("bottomright",paste(col_legend[,1]," : ",col_legend[,2],sep=""),pch=1,col = 1:dim(col_legend)[1])
					png_bool <- png_bool +1
					last_plot <- last_plot +1

					if(png_bool == 4){
						dev.off()
						cat(paste("see : ",outputfile,nb_output,".png\n",sep=""))
						png_bool <- 0
						nb_page <- nb_page -1
						nb_output <- nb_output +1
					}
				}
				else{
					if(png_bool == 0){png(paste(outputfile,nb_output,".png",sep=""),  width=860, height=860); par(mfrow=c(2,2))}

					# creation du plot
					plot(tmp$compA, tmp$compB, main="plot components", xlab=col_names[indA], ylab=col_names[indB],col=tmp$col)
					legend("bottomright",paste(col_legend[,1]," : ",col_legend[,2],sep=""),pch=1,col = 1:dim(col_legend)[1])
					png_bool <- png_bool +1
					last_plot <- last_plot +1

					if(last_plot == nb_plot){
						dev.off()
						cat(paste("see : ",outputfile,nb_output,".png\n",sep=""))
						png_bool <- 0
						nb_page <- nb_page -1
					}
				}
			}
		}
	}
}

#####

# creation de plots montrant l'association d'individus aux composantes en fonction de l'etude desquels ils sont tires
plot.indiv.components.orig <- function(nb_comp=NULL, all_famfile=NULL, casesctrls_famfile=NULL, ctrls_famfiles=NULL, phenofile="pheno_strat.txt",outputfile="qqplot_components_orig"){
	
	if(!is.null(nb_comp) & !is.null(all_famfile) & !is.null(ctrls_famfiles)){

		# chargement des donnees
		all <- read.table(file=all_famfile)
#		for (i in 1:(length(ctrls_famfiles))){
#			tmp <- read.table(file=ctrls_famfiles[i])
#			if(i==1){ all$orig <- ifelse(all$V2 %in% tmp$V2,i+2,NA) }
#			else{ all$orig <- ifelse(all$V2 %in% tmp$V2,i+2,all$orig)}
#		}

		if (!is.null(casesctrls_famfile)){
			study <- read.table(file=casesctrls_famfile)
			all$orig <- ifelse(((all$V2 %in% study$V2) & (all$V6==1)),i+1,all$orig) 
			all$orig <- ifelse(((all$V2 %in% study$V2) & (all$V6==2)),i+2,all$orig) 
		}
		for (i in 1:(length(ctrls_famfiles))){
			tmp <- read.table(file=ctrls_famfiles[i])
			if(i==1){ all$orig <- ifelse(all$V2 %in% tmp$V2,i,NA) }
			else{ all$orig <- ifelse(all$V2 %in% tmp$V2,i,all$orig)}
		}
		strat <- read.table(file=phenofile, h=T)
		col_names <- colnames(strat[-c(1,2)])
		strat <- merge(strat, all, by.x="IID", by.y="V2")
		strat <- strat[-c(1,2)] 
		colnames(strat) <- NULL
		
		nb_plot <- ((nb_comp)*(nb_comp -1))/2
		nb_page <- ifelse(nb_plot%%4 == 0, nb_plot%/%4, (nb_plot%/%4)+1)
		nb_output <- 1 # compte le nombre de fichiers en sortie
		png_bool <- 0 # permet la gestion du nombre de fichiers png (en fonction du nombre de plots dejà mis sur le png)
		last_plot <-0

		for(indA in 1:(nb_comp -1) ){
			for(indB in (indA+1):nb_comp){
				compA <- strat[c(indA)]
				compB <- strat[c(indB)]
				col <- strat[c(length(strat))]
				head(col)
				# creation des parametres de la legende 
				if (!is.null(casesctrls_famfile)){
					names <- c("controls","cases",ctrls_famfiles)
					col_legend<- cbind(c(1:(length(ctrls_famfiles)+2)),names)
				}
				else{
					names <- c(ctrls_famfiles)
					col_legend<- cbind(c(1:(length(ctrls_famfiles))),names)
				}

				# creation de la matrice de donnees servant au plot
				tmp <- cbind(compA,compB,col)	
				if(nb_page > 1){
					if(png_bool == 0){png(paste(outputfile,nb_output,".png",sep=""), width=860, height=860); par(mfrow=c(2,2))}#on ouvre un nouveau png tous les 4 plots
					# creation du plot
					plot(tmp$compA, tmp$compB, main="plot components", xlab=col_names[indA], ylab=col_names[indB],col=tmp$col)
					legend("bottomright",paste(col_legend[,1]," : ",col_legend[,2],sep=""),pch=1,col = 1:dim(col_legend)[1])
					png_bool <- png_bool +1
					last_plot <- last_plot +1

					if(png_bool == 4){
						dev.off()
						cat(paste("see : ",outputfile,nb_output,".png\n",sep=""))
						png_bool <- 0
						nb_page <- nb_page -1
						nb_output <- nb_output +1
					}
				}
				else{
					if(png_bool == 0){png(paste(outputfile,nb_output,".png",sep=""),  width=860, height=860); par(mfrow=c(2,2))}

					# creation du plot
					plot(tmp$compA, tmp$compB, main="plot components", xlab=col_names[indA], ylab=col_names[indB],col=tmp$col)
					legend("bottomright",paste(col_legend[,1]," : ",col_legend[,2],sep=""),pch=1, col = 1:dim(col_legend)[1])
					png_bool <- png_bool +1
					last_plot <- last_plot +1

					if(last_plot == nb_plot){
						dev.off()
						cat(paste("see : ",outputfile,nb_output,".png\n",sep=""))
						png_bool <- 0
						nb_page <- nb_page -1
					}
				}
			}
		}
	}
}

#####

# creation d'un unique plot avec une ligne horizontale ou verticale
plot.indiv.components.abline <- function(comp_number=NULL,famfile=NULL,phenofile=NULL,color="status",Hline=NULL,Vline=NULL){#comp_number=c(1,2) : component1 VS component2, Hline : ligne horizontale, Vline : ligne verticale
	if(!is.null(comp_number) & !is.null(famfile) & !is.null(phenofile)){
		# chargement des donnees
		fam <- read.table(file=famfile)
		strat <- read.table(file=phenofile, h=T)
		col_names <- colnames(strat[-c(1,2)])
		strat <- merge(strat, fam, by.x="IID", by.y="V2")
		strat <- strat[-c(1,2)] 
		colnames(strat) <- NULL

		if(color=="gender"){
			col <- strat[c(length(strat)-1)]
			# creation des parametres de la legende 
			col_legend<- cbind(c(1,2),c("male","female"))
		}
		if(color=="status"){
			col <- strat[c(length(strat))]
			col_legend<- cbind(c(1,2),c("controls","cases"))
		}
		# creation du plot
		compA <- strat[comp_number[1]]
		compB <- strat[comp_number[2]]
		tmp <- cbind(compA,compB,col)
		plot(tmp$compA, tmp$compB, main="plot components", xlab=paste("C",comp_number[1],sep=""), ylab=paste("C",comp_number[2],sep=""),col=tmp$col)
		legend("bottomright",paste(col_legend[,1]," : ",col_legend[,2],sep=""),pch=1,col = 1:dim(col_legend)[1])
		abline(h=Hline, v=Vline,col=4)
	}
}

#####

# creation d'un tableau bilan des donnees d'association aux composantes
summary.qassoc <- function(nb_comp=NULL, assocfile="pheno", outputfile="summary-qassoc.txt"){

	if(!is.null(nb_comp)){
		# chargement des donnees
		for( i in 1:nb_comp){

			Ctmp <- read.table(file=paste(assocfile,".C",i,".qassoc",sep=""),h=T)

			if(i == 1){ 
				table <- Ctmp[c(1,2,3,9,8)]
				names(table)[4] <- "PC1"
				names(table)[5] <- "TC1"
				table$Tsq1 <- (table$TC1)^2
				names(table)[6] <- "Tsq1"
				ind <- 6
			}
			else{
				table <- cbind(table,Ctmp[c(9,8)])
				ind <- ind+2

				# ajout des statistiques au carre
				Tsq <-(table[c(ind)])^2
				table <- cbind(table,Tsq)
			
				# on nomme les colonnes
				T_name <- paste("TC",i,sep="")
				P_name <- paste("PC",i,sep="")
				Tsq_name <- paste("Tsq",i,sep="")
				names(table)[ind-1] <- P_name
				names(table)[ind] <- T_name
				names(table)[ind+1] <- Tsq_name
				ind <- ind+1
			}
		}

		# sauvegarde des donnees
		write.table(table, file=outputfile, row.names=F,quote=F)
		cat(" see :",outputfile," \n")
	}
	else { cat("missing nb_comp parameter")}
}

#####

# creation de qqplot d'association aux composantes
qqplot.components<- function(nb_comp=NULL, IC=NULL, summaryfile="summary-qassoc.txt", outputfile_sq="qqplot_statsq_assoc_components", outputfile_p="qqplot_pval_assoc_components"){#outputfile_sq = nom du fichier des qqplots avec la stat au carre, outputfile_p = nom du fichier des qqplots des p-values
	
	#source("qqmultiple_function.R")
	data <- read.table(file=summaryfile ,h=T)
	data <- na.omit(data)
	library(car)
 
	nb_page <- ifelse(nb_comp%%4 == 0, nb_comp%/%4, (nb_comp%/%4)+1)	
	Pcol <- 4
	Tsqcol <- 6

	for(i in 1:nb_page){

		# plots pour la statistique au carre
		png(paste(outputfile_sq,i, ".png", sep=""),width=860, heigh=860)
		par(mfrow=c(2,2))
		j <- 1

		if(nb_page > 1){
			if(i < nb_page){
				for(j in 1:4){
					qqPlot(data[c(Tsqcol)], distribution="chisq", df=1, main="qqplot squared statistic", xlab="Expected values", ylab="Observed test statistic", grid=F)
					legend("bottomright",names(data)[Tsqcol],pch=1,col = 1)
					Tsqcol <- Tsqcol+3
				}
			}
			else{
				for(k in 1:(nb_comp-(4*(nb_page-1)))){
					qqPlot(data[c(Tsqcol)], distribution="chisq", df=1, main="qqplot squared statistic", xlab="Expected values", ylab="Observed test statistic", grid=F)
					legend("bottomright",names(data)[Tsqcol],pch=1,col = 1)
					Tsqcol <- Tsqcol +3
				}
			}
		}
		else{
			for(y in 1:nb_comp){
				qqPlot(data[c(Tsqcol)], distribution="chisq", df=1, main="qqplot squared statistic", xlab="Expected values", ylab="Observed test statistic", grid=F)
				legend("bottomright",names(data)[Tsqcol],pch=1,col = 1)
				Tsqcol <- Tsqcol+3
			}
		}
		dev.off()	
		cat(" > > > see : ",paste(outputfile_sq,i, ".png", sep=""),"\n")

		# plots pour les p-values
		data[data == 0] <- NA
		data <- na.omit(data)

		png(paste(outputfile_p, i, ".png", sep=""),width=860, heigh=860)
		par(mfrow=c(2,2))
			
		if(nb_page > 1){
			if(i < nb_page){
				for(j in 1:4){
					Pname <- names(data)[Pcol]
					plot.QQs(data1=na.omit(data), pobs.name=Pname, IC=IC)
					Pcol <- Pcol + 3
				}
			}
			else{
				for(k in 1:(nb_comp-(4*(nb_page-1)))){
					Pname <- names(data)[Pcol]
					plot.QQs(data1=na.omit(data), pobs.name=Pname, IC=IC)
					Pcol <- Pcol + 3
				}
			}
		}
		else{
			for(y in 1:nb_comp){
				Pname <- names(data)[Pcol]
				plot.QQs(data1=na.omit(data), pobs.name=Pname, IC=IC)
				Pcol <- Pcol + 3
			}
		}
		dev.off()
		cat(" > > > see : ",paste(outputfile_p, i, ".png", sep=""),"\n")
	}
	
}


######################### ANALYSE PAR PLAQUE #############################
###(creation des qqplots de comparaison de plaque les unes contre les autres et manhattan plot par plaque)###

# creation des fichiers cluster - a faire pour les cas et les controles
create_cluster_file <- function(inbfile,plate_file, num_Pdel=NULL, status) # num_Pdel : numero de plaque à supprimer, status : 1 for controls, 2 for cases
{
	# chargement et preparation des donnees
	plates <- read.table(file=plate_file,h=T)
	names(plates)[2] <- "ADN"
	names(plates)[3] <- "NPLAK"
	if(status == 2){
		system(paste("plink --noweb --bfile ",inbfile," --filter-cases --make-bed --out cases",sep=""))
		info <- read.table(file="cases.fam")
	}
	if(status == 1){
		system(paste("plink --noweb --bfile  ",inbfile," --filter-controls --make-bed --out controls",sep=""))
		info <- read.table(file="controls.fam")
	}

	info <- info[-c(3,4,5,6)]
	names(info)[1] <- "FID"
	names(info)[2] <- "ADN"

	# fusion des donnees
	cluster <- merge(info,plates, by=c("ADN"))
	cluster <- cluster[-c(3)]
	names(cluster)[1] <- "IID"
	cluster<- cluster[c(2,1,3)]
	names(cluster)[3] <- "NPLAK"	

	if(!is.null(num_Pdel)){
		if(length(num_Pdel) != 0){
			# creation du fichier cluster 
			for(i in 1:length(num_Pdel)){
				tmp_cluster <- cluster[cluster$NPLAK != num_Pdel[i],]
				cluster <- tmp_cluster
			}
		}
	}

	# enregistrement des donnees
	if(status == 2) {
	write.table(cluster,file="cluster-cases.dat",quote=F,row.names=F)
	cat("see results : cluster-cases.dat\n\n")
	}
	if(status == 1){
	write.table(cluster,file="cluster-ctrls.dat",quote=F, row.names=F)
	cat("see results : cluster-ctrls.dat\n\n")
	}
}

#####

# creation du fichier de phenotype - a faire pour les cas et les controles
create_pheno_file <- function(cluster_file, status) #status : 1 for controls, 2 for cases
{
	# chargement des donnees
	cluster <- read.table(file=cluster_file,h=T)
	plate_list <-sort(unique(cluster$NPLAK))
	x <- length(cluster$FID)
	y <- length(plate_list)

	# creation de la matrice des donnees
	mat <- matrix(1,nrow=x,ncol=y)
	nplak <- cluster$NPLAK

	# remplissage de la matrice de donnees
	for (i in 1:x){
	j=1
	repeat{
	if(nplak[i]==plate_list[j])break else{j=j+1;}
	}
	mat[i,j] <- 2
	}

	# mise en forme de la matrice
	pheno <- cbind(cluster[-c(3)],mat)

	# nommer les colonnes
	names <- colnames(cluster[-c(3)])
	col_names <- c(names,plate_list)
	colnames(pheno) <- col_names
	num <- colnames(pheno[-c(1,2)])

	# enregistrement des donnees
	if(status == 2) {
	write.table(pheno,file="pheno_cases.txt",quote=F, row.names=F, col.names=T)
	cat("see results : pheno_cases.txt\n")
	write.table(num,file="header-pheno-file1.txt")
	}
	if(status == 1){
	write.table(pheno,file="pheno_ctrls.txt",quote=F, row.names=F, col.names=T)
	cat("see results : pheno_ctrls.txt\n")
	write.table(num,file="header-pheno-file2.txt")
	}
}

#####

# creation des fichiers d'association sous plink
association <- function(status=NULL)
{
	if(status == 1){
		# pour les contrôles
		system(paste("plink --noweb --bfile controls --pheno pheno_ctrls.txt --all-pheno --assoc --out ctrls_plate",sep=""))
		cat("create controls.assoc \n")
	}

	if(status == 2){
		# pour les cas
 		system(paste("plink --noweb --bfile cases --pheno pheno_cases.txt --all-pheno --assoc --out cases_plate",sep=""))
		cat("create cases.assoc \n")
	}
}

#####

# sert pour la fonction info_chisq_pval (renvoie le nom de la colonne qui correspond au minimum)
name <- function(data){
	ind <- which.min(data)
	val <- names(data)[ind]
}

# sert pour la fonction info_chisq_pval (calcul du lambda median)
lambda.median <- function(data){
	res <- (median(data))/0.456
}

# creation d'un tableau d'informations sur les chi² et les p-values
info_chisq_pval <- function(p.names="P", status=NULL){ #status : 1 for controls, 2 for cases
	
	if(status == 2){
		# chargement du fichier de phenotype
		assoc_name <- "cases_plate"
		num <- read.table(file="header-pheno-file1.txt")
	}
	if(status == 1){
		# chargement du fichier de phenotype
		assoc_name <- "ctrls_plate"
		num <- read.table(file="header-pheno-file2.txt")
	}	

	# chargement des donnees de plaque
	for ( i in 1:length(p.names)) {
		name <- paste(assoc_name, num[i,1], sep=".")
		name <- read.table(file=paste(name,"assoc", sep="."),h=T)
		if(i == 1) { 
			pval <- name[c(1,2,3)]
			chisq <- name[c(1,2,3)]
		}
		tmp_pval <- name$P
		tmp_chisq <- name$CHISQ

		#  fusion des donnees
		pval <- cbind(pval,tmp_pval)
		chisq <- cbind(chisq,tmp_chisq)
	}

	#  on nomme les colonnes
	col_names <- c("CHR","SNP","BP",p.names)
	colnames(pval) <- col_names
	colnames(chisq) <- col_names
	
	pval <- na.omit(pval)
	pval_complete <- pval
	chisq_complete <- chisq

	# moyenne des chisq par plaque 
	chisq_tmp <- chisq[-c(1,2,3)]
	mean_chisq <- unlist(apply(X=chisq_tmp,na.rm=T , MARGIN=2, FUN=mean))

	# mediane des chisq par plaque
	lambda_median_chisq <- unlist(apply(X=chisq_tmp , MARGIN=2, FUN=lambda.median))

	# minimum des p-value par plaque
	pval_tmp <- pval[-c(1,2,3)]
	min_pval <- unlist(apply(X=pval_tmp,na.rm=T , MARGIN=2, FUN=min))

	# creation d'un fichier info_chisq_pval
	info_chisq_pval <- cbind(p.names,mean_chisq,lambda_median_chisq,min_pval)	

	# ajout de la p-value minimum et du nom de la plaque associee	
	pval_tmp <- pval[-c(1,2,3)]
	pval$P_min <- apply(X=pval_tmp, MARGIN=1, FUN=min)
	pval$P_name <- apply(X=pval_tmp, MARGIN=1, FUN=name)

	# enregistrement des donnees
	if(status == 2){
		write.table(pval_complete,file="p-values-complete-cases.txt",quote=F,row.names=F)
		write.table(chisq_complete,file="chisq-complete-cases.txt",quote=F,row.names=F)

		write.table(pval,file="p-values-cases.txt", quote=F, row.names=F)
		cat("see p-values-cases.txt\n")
		write.table(info_chisq_pval, file="info_chisq_pval_cases.txt", quote=F, row.names=F,sep="\t")
		cat("see info_chisq_pval_cases.txt\n")
	}
	if(status ==1){
		write.table(pval_complete,file="p-values-complete-controls.txt",quote=F,row.names=F)
		write.table(chisq_complete,file="chisq-complete-controls.txt",quote=F,row.names=F)

		write.table(pval,file="p-values-controls.txt", quote=F, row.names=F)
		cat("see p-values-controls.txt\n")
		write.table(info_chisq_pval, file="info_chisq_pval_controls.txt", quote=F, row.names=F,sep="\t")
		cat("see info_chisq_pval_controls.txt\n")
	}
}

#####

# creation d'un QQplot par plaque - a faire pour les cas et pour les controles
qqplot.p<- function (p.names="P", nb_by_graph=4, IC=NULL, status=NULL, outputfile=NULL) #nb_by_graph : nb de plaque par graphique ; status : 1 for controls, 2 for cases
{
	if((nb_by_graph > length(p.names)) | (nb_by_graph==0)) {
		stop("ERROR : incorrect nb_by_graph \n")
	}
	else{

		#  visualisation des p-values des plaques 
		nb_plate <- length(p.names)
		nb_plot <- ifelse(nb_plate%%nb_by_graph == 0, nb_plate%/%nb_by_graph, (nb_plate%/%nb_by_graph)+1)
		nb_page <- ifelse(nb_plot%%4 == 0, nb_plot%/%4, (nb_plot%/%4)+1)
		x<-1
		y<-nb_by_graph

		if(status == 2){
			if(is.null(outputfile)){outputfile <- "qqplot-plaque-cases"}
			# chargement des donnees
			pval <- read.table(file="p-values-cases.txt",h=T)
		}
		if(status == 1){
		if(is.null(outputfile)){outputfile <- "qqplot-plaque-controls"}
			# chargement des donnees
			pval <- read.table(file="p-values-controls.txt",h=T)
		}

		for (k in 1:nb_page){
				png(paste(outputfile,k,".png", sep=""), width=940, height=940, res=72)
				par(mfrow=c(2,2))		

				if(nb_page >1){
					if (k < nb_page){
						for (j in 1:4){ 
							plot.QQs(data1=pval,pobs.name=p.names[x:y],IC=IC)
							nb_plot <- nb_plot-1
							x<-x+nb_by_graph
							y <- y+nb_by_graph
						}
					}
					else{
						if(y > length(p.names)){y<-length(p.names)}	
						for (i in 1:(nb_plot)){
							plot.QQs(data1=pval,pobs.name=p.names[x:y],IC=IC)
							nb_plot <- nb_plot-1
							x <- x+nb_by_graph
							if((y+nb_by_graph) < length(p.names)){y <- y+nb_by_graph}
							else if(x <= length(p.names)){
								y <- length(p.names)
								plot.QQs(data1=pval,pobs.name=p.names[x:y],IC=IC)
								nb_plot <- nb_plot-1
								break
							}
						}
					}
				}
				else{
					for (h in 1:nb_plot){ 
						plot.QQs(data1=pval,pobs.name=p.names[x:y],IC=IC)
						x <- x+nb_by_graph
						if((y+nb_by_graph) < length(p.names)){y <- y+nb_by_graph}
						else if(x <= length(p.names)){
							y<- length(p.names)
							plot.QQs(data1=pval,pobs.name=p.names[x:y],IC=IC)
							break
						}
					}
				}
				dev.off()
		}
	}
}

#####

# creation d'un manhattan plot pour chaque plaque
mhtplot.p <- function (p.name=NULL, status=NULL,title=NULL,outputfile=NULL) #status : 1 for controls, 2 for cases
{
	# chargement du package "gap"
	library(gap)

	# gestion des erreurs de saisie du paramètre p.name
	if((is.null(p.name)) || (length(p.name) >1)){
		print("error : p.name needs one plak name")
	}
	else{
		if(status == 2){
			pval <- read.table(file="p-values-cases.txt", h=T)
			chr_pos <- pval[c(1,3)]
		}
		if(status == 1){
			pval <- read.table(file="p-values-controls.txt", h=T)
			chr_pos <- pval[c(1,3)]
		}

		# mise en page du mhtplot
		#colors <- c("red","blue","green","cyan","magenta","gray","red","blue","green","cyan","magenta","gray","red","blue","green","cyan","magenta","gray","red","blue","green","cyan","magenta","gray")
     		colors <- rep(c("blue","green"),12)

		# recuperation des donnees pour le manhattan plot
		p <- with(pval, eval(parse(text = p.name)))# recuperation de la colonne correspondant a la plaque d'interet
		pobs <- cbind(chr_pos,p)

		# creation du manhattan plot
		if(is.null(outputfile)) {outputfile <- paste("mhtplot",p.name,sep="_")}
		png(paste(outputfile,".png",sep=""),res=55, width=1200,height=650)

		mhtplot(data=pobs,control=mht.control(usepos=TRUE,,colors=colors,cutoffs=c((-log10(5e-08)),6)),pch=19)
		axis(2,pos=2)
		par(cex=1.5, cex.axis=1.5) 
		if(!is.null(title)) {title(title)}

		dev.off()
	}
}

#####

# mise à zéro de certains marqueurs jugés trop différents au niveau des fréquences alléliques
obligatory_missing_snp <- function(inbfile=NULL,pvalFile=NULL,th=NULL,status=NULL,method=NULL,platefile=NULL){
	if(!is.null(inbfile) & !is.null(pvalFile) & !is.null(th) & !is.null(status) & !is.null(method)){
		#chargement des donnees
		pval <- read.table(file=pvalFile,h=T)

		zerofile <- NULL
		clusterfile <- NULL

		# creation d'une matrice dont tous les marqueurs superieurs au seuil sont mis comme "NA"
		data <-  apply(pval[4:(dim(pval)[2]-2)],2,function(data) ifelse(as.numeric(data) < th,data,"NA"))
		#data <-  apply(pval[-c(1,2,3)],2,function(data) ifelse(as.numeric(data) < th,data,"NA"))
		# ajout des "rsID" et creation d'un objet de type "data frame"
		data2 <- as.data.frame(cbind(as.character(pval$SNP),data))
		col_names <- colnames(data2)

		if(method == 1){

			# prepartation des fichiers cluster ".clst" (format plink)
			if(status == 1){			
				header <- read.table(file="header-pheno-file2.txt",h=T)
				cluster <- read.table(file="cluster-ctrls.dat",h=T)
			}
			if(status == 2){
				header <- read.table(file="header-pheno-file1.txt",h=T)
				cluster <- read.table(file="cluster-cases.dat",h=T)
			}
			corresp <- cbind(col_names[2:length(col_names)],header) # correspondance nom plaque / numero plaque
		}else if(method == 2 & !is.null(platefile)){
			cluster <- read.table(file=platefile)
			platenames <- read.table(file="plates-names.txt")
		}

		for( i in 2:dim(data2)[2]){ # pour chaque plaque

			# creation des fichiers ".zero" (format plink) : liste de marqueurs à mettre a zero pour chaque plaque + phenotype(plaque1...)
			tmpdata2 <- data2[c(1,i)] # rsIBD plaque
			list_to_del <- tmpdata2[tmpdata2[c(2)] != "NA","V1"]# liste des SNPs dont la p-valeur <seuil 
#			col_names <- colnames(data2)
			if(length(list_to_del) !=0){ # il y a des SNPs en dessous du seuil

				phenotype <- rep(as.character(platenames[i-1,1]),length(list_to_del))# creation de la colonne phenotype du fichier ".zero"
				tmpzerofile <- as.data.frame(cbind(as.character(list_to_del),phenotype)) # fichier ".zero"

				# creation des fichiers cluster ".clst" (format plink)
				if(method==1){
					tmpclusterfile1 <- cluster[cluster[c(3)] == corresp[i-1,2] , c("FID","IID")] # recuperation des individus correspondant a la plaque
				} else if(method==2){
					phenotype_substr <- substr(phenotype[1],2,nchar(phenotype[1])) # on enlève le premier caractère
					tmpclusterfile1 <- cluster[cluster[c(3)] == phenotype[1] | cluster[c(3)] == phenotype_substr, c("V1","V2")] 
				}
				phenotype <- rep(platenames[i-1,1],dim(tmpclusterfile1)[1])# creation de la colonne phenotype du fichier ".clst"
				tmpclusterfile2 <- as.data.frame(cbind(tmpclusterfile1,phenotype)) # fichier ".clst"				

				# concatenation des donnees 
				zerofile <- rbind(zerofile,tmpzerofile)
				clusterfile <- rbind(clusterfile,tmpclusterfile2)
 				
			}
		}
					# enregistrement des donnees
			write.table(zerofile,file=paste("snps_under_",th,".zero",sep=""),quote=F,row.names=F,col.names=F)
			write.table(clusterfile,file=paste("snps_under_",th,".clst",sep=""),quote=F,row.names=F,col.names=F)

			#  mise a zero des SNPs dont la p-valeur < seuil
			system(paste("plink --noweb --bfile ",inbfile," --oblig-missing snps_under_",th,".zero --oblig-clusters snps_under_",th,".clst --make-bed --out oblig_missing",sep=""))
	}
}




######################## ANALYSE D'ASSOCIATION ############################
# creation des plots de résultats de l'analyse d'association
plots.association.analysis <- function(associationfile=NULL, IC=NULL, qqplot_title="QQPlot", mhtplot_title="Manhattan plot", out_qqplotfile="qqplot", out_mhtplotfile="mhtplot")
{
	# mise en forme des données pour le mhtplot et le qqplot 
	if(!is.null(associationfile) && !is.null(IC)){
		logistic <- read.table(file=associationfile,h=T)
		logisticADD <- logistic[logistic$TEST=="ADD",]
	}
	else{
		cat("ERROR : associationfile or IC not found")
	}
	# qqplot final
	png(paste(out_qqplotfile,".png",sep=""))
	plot.QQs(data1=logisticADD,pobs.name="P",IC=IC)
	title(qqplot_title) 
	dev.off()
	cat("see",out_qqplotfile,"\n")

	# mhtplot final
	library(gap) 
	pheno <- logisticADD[c("CHR","BP","P")] 
	colors <- rep(c("blue","green"),12)
	png(paste(out_mhtplotfile,".png",sep=""),res=55, width=1200,height=650)
	mhtplot(data=pheno, control=mht.control(usepos=TRUE,colors=colors,cutoffs=c((-log10(5e-08)),6)), pch=19)
	axis(2,pos=2) 
	par(cex=1.5, cex.axis=1.5) 
	title(mhtplot_title)
	dev.off()
	cat("see",out_mhtplotfile,"\n")
}


############################# FONCTION qqmultiple ##########################
plot.QQs <- function(	data1 = dat, 
			subset = NULL, #"ma_name == \"Hind\" & ma_type == \"affy\"",
			pobs.name = "pvalue", #c("pvalTest1", "pvalTest2")
			lim.X = NULL, ##c(xmin, xmax)
			lim.Y = NULL, ## c(ymin, ymax), 
			titre = NULL,
			couleur = T,	
			IC=NULL){

qqplot <- function (x, y, plot.it = TRUE, xlab = deparse(substitute(x)), 
    ylab = deparse(substitute(y)), add.to.plot = FALSE, ...) 
{
    sx <- sort(x)
    sy <- sort(y)
    lenx <- length(sx)
    leny <- length(sy)
    if (leny < lenx) 
        sx <- approx(1:lenx, sx, n = leny)$y
    if (leny > lenx) 
        sy <- approx(1:leny, sy, n = lenx)$y
    if (plot.it)
      if(add.to.plot)
        points(sx, sy, ...)
      else
        plot(sx, sy, xlab = xlab, ylab = ylab, ...)
    invisible(list(x = sx, y = sy))
}

confidence_intervals<- function (x, y)
{
    #x vector of pval
	#y confidence interval
	
	observed <- -(log10(x))
	N <- length(x) # number of p-values
	uniform <- -log10(1:N/(N+1))
	MAX <- max(c(observed, uniform))
	


	##Calcul d'interval de confiance
	
	q05 <- rep(0,N)
	q95 <- rep(0,N)
	
	for(i in 1:N)
	{
		#La k-ème statistique d'ordre d'un n-échantillon de lois uniformes U(0,1] suit la loi Beta(k,n-k+1) 
		q05[i] <- qbeta(1-y,i,N-i+1)
		q95[i] <- qbeta(y,i,N-i+1)
		
	}

	points(uniform,-log10(q95),col="black",type="l", lty=2, lwd=2)
	points(uniform,-log10(q05),col="black",type="l", lty=2, lwd=2)
	
}

print(dim(data1))
											
	if(!is.null(subset)) data1 <- data1[with(data1, eval(parse(text = subset))), ]
	
print(dim(data1))

numPval <- vector(length= length(pobs.name))

	#scale
	if(-log10(min(data1[,which(names(data1) %in% pobs.name)],na.rm=T)) == -log10(0)) min <- -log10(1e-200)
		else min <-  -log10(min(data1[,which(names(data1) %in% pobs.name)],na.rm=T))

#	if(-log10(max(data1[,which(names(data1) %in% pobs.name)],na.rm=T))== -log10(1)) max <- 1
#		else max <-  -log10(max(data1[,which(names(data1) %in% pobs.name)],na.rm=T))

#	if(is.null(lim.Y))  testY <- c(-log10(max(data1[,which(names(data1) %in% pobs.name)],na.rm=T)), -log10(min(data1[,which(names(data1) %in% pobs.name)],na.rm=T)))
	if(is.null(lim.Y))  testY <- c(-log10(max(data1[,which(names(data1) %in% pobs.name)],na.rm=T)), min)
		else testY <- lim.Y

	if(is.null(lim.X)) testX <- testY
		else testX <- lim.X

for(i in 1:(length(pobs.name))){

 cat(paste("\n**********\ngraphe", i, ":", pobs.name[i], "\n**********\n"))

	pobs <- with(data1, eval(parse(text = pobs.name[i])))

	if(is.factor(pobs)) pobs <- as.numeric(levels(pobs))[pobs]

 cat(paste("nb p-values renseignees :", sum(!is.na(pobs))))

	numPval[i] <- sum(!is.na(pobs))

	#informations
	if(is.null(lim.X))	scaleX <- c(-log10(max(pobs, na.rm = T)), -log10(min(pobs, na.rm = T)))
		else scaleX <- lim.X
	
	if(is.null(lim.Y))  scaleY <- c(-log10(max(pobs, na.rm = T)), -log10(min(pobs, na.rm = T)))
		else scaleY <- lim.Y
	
cat("\nlimites X : "); print(scaleX)
cat("limites Y : "); print(scaleY)
	
	pobsmin <- min(pobs, na.rm = T)
	pobsmax <- max(pobs, na.rm = T)

cat(paste("pobsmin :", pobsmin))
cat(paste("\npobsmax :", pobsmax))	

cat(paste("\n-log10(pobsmin) :", -log10(pobsmin)))	
cat(paste("\n-log10(pobsmax) :", -log10(pobsmax), "\n"))

	qqplot( -log10((1:length(pobs)/(length(pobs) + 1))), #rang

		-log10(pobs), 
		main = paste("QQ-plot", subset, titre, sep = "\n"), 
		ylab="-log10(observed p-value)", 
		xlab="-log10(uniform p-value)", 
		xlim = testX,
  		ylim = testY , 
		pch = i,
		cex = 0.6,
		col = ifelse(couleur, i+1 ,1),
		add.to.plot = ifelse(i > 1, T, F)

)

	abline(0, 1, col="black", lwd=2)

}

# confidence interval
if (!is.null(IC))
{
		confidence_intervals( pobs, y=IC )  
}	

#legend
legend("bottomright", paste(pobs.name, " (N = ",  numPval, ")", sep = ""), pch = 1:length(pobs.name), col = 2:(length(pobs.name)+1))


}
